Learning Outcomes

• These lectures will teach you how to:
  • Identify when a random variable follows a binomial distribution, and what its parameters are
  • Understand what a normal distribution is, and what its parameters are
  • Compute important numeric features of these distributions using the d<dist>, p<dist>, and q<dist> commands

Preliminaries

• You will need to install the egg package for a couple isolated examples in these lecture notes, with install.packages("egg").

1. Download the file week07-distributions.Rmd into STAT240/lecture/week07-distributions.

Assumptions

• In this class, we will primarily derive distributions the second way, by making reasonable assumptions.

• For example, forgetting all of the numbers involved, let’s try to work out the probability distribution of how many red lights I will encounter out of five total lights on the way to school.

• In reality, the traffic lights are often coordinated with each other, creating a complex pattern of dependency with what lights I make it through and which I stop at.

• This makes it very hard, if not impossible, to determine the “real” probability distribution.

• However, if I make two reasonable assumptions, it becomes very easy to determine the distribution.

  • Assumption #1: The state of each light as I approach it does not affect and is not affected by any other lights - they are independent.
  • Assumption #2: All lights have the same probability p of being green as I approach.

These assumptions are definitely simplifications of reality, as we just stated lights depend on each other in a complex way. However, they are reasonable simplifications.

• By making these assumptions, I can determine the whole probability distribution.

• For example, the probability of hitting five straight greens is \(p*p*p*p*p = p^5\), the probability of hitting five straight reds is \((1-p)^5\), and the probability of hitting four greens and a red is \(5 * p^4 * (1-p)\) (more on these formulas later)!

• The calculations above require the assumptions to be true.

• The assumptions, in the most technical sense, are not true. Many assumptions will never be “technically” true, and that’s okay. If we are willing to accept the assumptions as “true enough”, then we can actually start making some useful progress.

The BINS Assumptions

To assess whether a random variable follows a binomial distribution, check if it meets the four BINS assumptions:

• B = Binary Outcomes: Each individual trial may be considered either a success or failure. (In particular, it is not returning a number other than 0 or 1.)

• I = Independent: Each individual trial does not affect the outcome of another.

• N = Fixed Sample Size: The number of trials is pre-specified, and does not depend on the outcome of those trials.

• S = Same Probability: Each individual trial has the same probability of success as any other.

All of these criteria must be met for a random variable to be binomial. In other words, if any fail, your variable is not binomial.

• We can verify that \(X\) from above, the number of heads meets these criteria:

  • B: Binary Outcomes: Each trial either succeeds (heads) or fails (not heads)
  • I: Independent: The result of any coin flip does not affect the other coin flips.
  • N: Fixed Sample Size: We pre-specified we were going to flip the coin 10 times.
  • S: Same probability: Each coin flip has probability 0.5 of turning up heads. The coin does not get “better” or “worse” at turning up heads.

• Another example of a binomial random variable, let’s call this \(D\), is the number of times I roll a 4 from five rolls of a fair six-sided die.

  • We would say that \(D \sim Binom(5, 1/6)\).

  • The last three assumptions; independence, fixed sample size, and same probability, are clearly met for the same reasons as the previous example.

  • However, it is tempting to say that this is not a binary outcome, because a single roll of the die can produce any of six outcomes, whereas the coin was heads or tails.

  • This reasoning is tempting, but forgets the definition of our random variable.

  • If all we care about is whether or not we got a 4, each trial is binary; we either got a 4 or we didn’t.

1. Outcome of the random process (five rolls of a fair six-sided die): 4-1-2-4-5

2. How you could view this as a binary outcome: Success - Fail - Fail - Success - Fail

3. Value of the random variable: \(D = 2\)

• Two examples of random variables which are not binomial:

• The raw measurement of a continuous variable, such as the height of a person, is not binomial.

  • Because a binomial variable must be the number of successes from a certain number of trials, it can only take on integer values. A binomial variable will never take on a value like 50.5.
  • You would have to make each numeric trial into a success-fail, such as measuring “the number of people out of 20 who are at least 5’5”“.

Background note: In basketball, a free throw is an un-defended shot from about fifteen feet that you can take standing still, as compensation if the other team fouls you.

• You take free throws until you make one. Let \(F\) be the number of free throws you take in total.
  • While it is true that each trial is make-or-miss (B), your free throws don’t affect each other (I), and each one has the same basic probability of going in (S), the N: Fixed Sample Size criteria is not met.
  • We are not prespecifying how many trials we are going to undertake; what would we use for \(n\) in \(F \sim Binom(n, p)\)?

The Binomial Probability Distribution Informal Proof

We will not cover this proof in lecture nor will we ever test you on it. It is only here for the curious, and because Cameron put all that effort into writing it and thought it was cool.

• Consider rolling a fair six-sided die three times, and let \(D\) be the number of 1’s you observe. We know \(D \sim Binom(3, 1/6)\).

• What is \(P(D = 1)\)? (In other words, what is the probability that three rolls of a fair six-sided die will yield EXACTLY one “1”?)

  • Not an easy thing to compute in your head!

• There are \(2^n = 2^3= 8\) possible outcomes, when each trial is success/fail:

# Let TRUE represent a success (1), FALSE represent a failure (Not a 1)
three_rolls = tibble(
  roll1 = c(F, T, F, F, T, T, F, T),
  roll2 = c(F, F, T, F, T, F, T, T),
  roll3 = c(F, F, F, T, F, T, T, T)
)
three_rolls

## # A tibble: 8 × 3
##   roll1 roll2 roll3
##   <lgl> <lgl> <lgl>
## 1 FALSE FALSE FALSE
## 2 TRUE  FALSE FALSE
## 3 FALSE TRUE  FALSE
## 4 FALSE FALSE TRUE 
## 5 TRUE  TRUE  FALSE
## 6 TRUE  FALSE TRUE 
## 7 FALSE TRUE  TRUE 
## 8 TRUE  TRUE  TRUE

• Exactly three of these outcomes result in \(D = 1\).

three_rolls = three_rolls %>% 
  mutate(D = roll1 + roll2 + roll3)

three_rolls

## # A tibble: 8 × 4
##   roll1 roll2 roll3     D
##   <lgl> <lgl> <lgl> <int>
## 1 FALSE FALSE FALSE     0
## 2 TRUE  FALSE FALSE     1
## 3 FALSE TRUE  FALSE     1
## 4 FALSE FALSE TRUE      1
## 5 TRUE  TRUE  FALSE     2
## 6 TRUE  FALSE TRUE      2
## 7 FALSE TRUE  TRUE      2
## 8 TRUE  TRUE  TRUE      3

• This does not mean that \(P(D=1) = 3/8\), because not every outcome has the same probability!

• Since the trials are independent, the probability of a given outcome occurring is the product of the probabilities of each trial.

• For example, Success - Fail - Fail in this case has probability of 1/6 * 5/6 * 5/6.

## # A tibble: 8 × 6
##   roll1 roll2 roll3 prob_str         prob     D
##   <lgl> <lgl> <lgl> <chr>           <dbl> <int>
## 1 FALSE FALSE FALSE 5/6 * 5/6 * 5/6 0.579     0
## 2 TRUE  FALSE FALSE 1/6 * 5/6 * 5/6 0.116     1
## 3 FALSE TRUE  FALSE 5/6 * 1/6 * 5/6 0.116     1
## 4 FALSE FALSE TRUE  5/6 * 5/6 * 1/6 0.116     1
## 5 TRUE  TRUE  FALSE 1/6 * 1/6 * 5/6 0.023     2
## 6 TRUE  FALSE TRUE  1/6 * 5/6 * 1/6 0.023     2
## 7 FALSE TRUE  TRUE  5/6 * 1/6 * 1/6 0.023     2
## 8 TRUE  TRUE  TRUE  1/6 * 1/6 * 1/6 0.005     3

• Observe that each case which has the same number of successes (the same value of \(D\), the same number of 1’s observed) has the same probability of happening!

• That’s because the order we multiply in doesn’t matter. For example, in the second and third row, \(\frac{1}{6}*\frac{5}{6}*\frac{5}{6} = \frac{5}{6}*\frac{1}{6}*\frac{5}{6}\).

• Both of those examples evaluate to \(\frac{1}{6}^1*\frac{5}{6}^2\). Notice that \(\frac{1}{6}\) is the probability of success of a single trial, and \(\frac{5}{6}\) is the “probability of failure”: \(1 - \frac{1}{6}\).

• Furthermore, the “1” exponent is the number of successes, and the “2” exponent is the number of failures (3 trials - 1 success).

• More generally,

\[ \text{Prob. of one specific outcome} = \text{Prob. Success}^{\text{Num. Success}} * \text{Prob. Failure}^{\text{Num. Failure}} \]

• Now, because every trial either succeeds or fails, we only need to know the probability and number of successes to determine the probability and number of failures.

• Let’s say we have probability of success \(p\), and we are interested in an outcome with \(k\) successes. (There may be more than one outcome with \(k\) successes - right now we care about a specific one.)

\[ \text{Prob. of one specific outcome} = p^k * (1-p)^{n-k} \]

• We can use \(k = 1\), along with the known \(n = 3\) and \(p = 1/6\), to determine the probability of each of the individual outcomes with exactly 1 success.

p = 1/6
k = 1
n = 3
round(p^k * (1-p)^(n-k), 3)

## [1] 0.116

# They match!
three_rolls %>% filter(D == 1)

## # A tibble: 3 × 6
##   roll1 roll2 roll3 prob_str         prob     D
##   <lgl> <lgl> <lgl> <chr>           <dbl> <int>
## 1 TRUE  FALSE FALSE 1/6 * 5/6 * 5/6 0.116     1
## 2 FALSE TRUE  FALSE 5/6 * 1/6 * 5/6 0.116     1
## 3 FALSE FALSE TRUE  5/6 * 5/6 * 1/6 0.116     1

• Finally, the probability of getting \(D = 1\) is the probability that we get ANY of these outcomes.

• Here, we have \(0.116 + 0.116 + 0.116 = 3 * 0.116 = 0.348\)!

• More generally, this is \(C * (p^k * (1-p)^{n-k})\), where \(C\) is the number of ways to get \(k\) successes from \(n\) trials, and the remainder of the term is the probability we just derived of getting an individual one of those outcomes.

• Here, \(C\) would be 3; there are three ways to get one success from three trials, namely, getting the success on your first, second, or third attempt.

• It is not immediately obvious how to obtain this number \(C\) mathematically.

• However, this number can be calculated more generally using something called the binomial coefficient.

• The binomial coefficient, \(\binom{n}{k}\), computes how many ways there are to get \(k\) successes from \(n\) trials.

\[ \binom{n}{k} = \frac{n!}{k!(n-k)!} \]

• The derivation of this coefficient is beyond the scope of this class, but a reasonable explanation for the curious can be found here.

Note: “!” represents the “factorial” expression, such that 3! = 3 times 2 times 1 = 6, 5! = 5 times 4 times 3 times 2 times 1 = 120, etc.

# The binomial coefficient can be calculated with choose(n, k)
choose(n, k)

## [1] 3

# Manually
factorial(n) / (factorial(k) * factorial(n - k))

## [1] 3

• Putting it all together, we get the crown jewel of the binomial distribution, the probability of getting \(k\) successes any way from \(n\) trials with success probability \(p\):

The Binomial Probability Distribution Result

\[ P(X = k) = \binom{n}{k} * p^k * (1-p)^{n-k} \ldots \text{ where } X \sim Binom(n, p) \]

… where \(n\) is the number of binomial trials, \(p\) is the (constant) probability of success of each trial, and \(k\) is the number of successes.

• You do not have to calculate this “manually” each time: R has a useful command to calculate this probability for you.

dbinom(k, size = n, prob = p), standing for “Density of the Binomial distribution”, will compute \(P(Binom(n, p) = k)\) for you.

• Let’s apply this formula to our above example; how likely is it that we will get a single “1” from three rolls of a fair, six-sided die?

n = 3
k = 1
p = 1/6

# From our theory:
choose(n, k) * p^k * (1-p)^(n-k)

## [1] 0.3472222

# This useful command calculates it for you! Definitely use this instead of the above code!
dbinom(k, n, p)

## [1] 0.3472222

• The probability of getting exactly one “1” from three rolls of a fair, six-sided die, written as \(P(Binom(3, 1/6) = 1)\), is roughly 34.7%.

• We can now enumerate the full probability distribution for \(D \sim Binom(3, 1/6)\):

# Reminder that these all should sum to 1!
n = 3
p = 1/6
tibble(
  x = 0:n,
  p = dbinom(x, n, p)
)

## # A tibble: 4 × 2
##       x       p
##   <int>   <dbl>
## 1     0 0.579  
## 2     1 0.347  
## 3     2 0.0694 
## 4     3 0.00463

Binomial Mean, Variance, and Standard Deviation

• Last week, we learned the general formula for the mean and variance of a probability distribution.

• The discrete formulas are shown below.

\[ E(X) = \sum x * P(X = x) \]

\[ Var(X) = \sum (x-\mu)^2*P(X=x) \]

• Those formulas will always work, but they are often inconvenient to compute. Last week’s examples were very small and friendly.

• Now that we know \(P(X = x) = \binom{n}{k} * p^k * (1-p)^{n-k}\) for any \(X \sim Binom(n, p)\), we can plug this into our mean and variance formulas.

• For example, to derive the mean of \(Binom(n, p)\), we can compute:

\[ E(X) = \sum_{k = 0}^n x * \binom{n}{k} * p^k * (1-p)^{n-k} \]

• No easy task, to be sure. Luckily for us, statisticians have already done this for us, and it simplifies to a very nice, intuitive result.

The mean of \(Binom(n, p)\) is n*p.

• For example, the mean of \(D \sim Binom(3, 1/6)\) is 3 * 1/6 = 0.5.

• Conducting a similar substitution and evaluation for the variance:

The variance of \(Binom(n, p)\) is \(\sigma^2 = np(1-p)\).

• The variance of \(D\) is \(3 * 1/6 * (1 - 1/6) \approx 0.417.\)

The standard deviation of \(Binom(n, p)\) is \(\sigma = \sqrt{np(1-p)}\).

• The standard deviation of \(D\) is \(\sqrt{Var(D)} = \sqrt{0.417} \approx 0.645\).

• These formulas are much simpler than the general forms!

• We will now investigate these formulas through a new concept called simulation.

# Try it out! Each value you get is the result of 10 coin flips.
rbinom(3, size = 10, prob = 0.5)

## [1] 7 7 4

random_data = rbinom(1000000, size = 100, prob = 0.5)
mean(random_data)

## [1] 49.99337

var(random_data)

## [1] 25.03594

dbinom(51, 100, 0.5)

## [1] 0.07802866

mean(random_data == 51)

## [1] 0.077844

# ggplot requires a dataframe, not a vector 
tibble(x = random_data) %>%  
  ggplot(aes(x)) +
  geom_histogram(
    aes(y = after_stat(density)), # Changes raw counts to observed probabilities
    center = 50, binwidth = 1, 
    color = "black", fill = "gray") +
  
  scale_y_continuous(expand = expansion(mult = c(0, 0.05)),
                     labels = scales::percent) +
  theme_classic() +
  
  labs(
    title = "Simulated Binomial Distribution",
    subtitle = "n = 100, p = 0.5",
    y = "Observed Probability"
  )

# Here's the real Binom(100, 0.5)!
gbinom(100, 0.5) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.05)),
                     labels = scales::percent) +
  xlim(c(25, 75)) +
  theme_classic() +
  labs(
    title = "Binom(100, 0.5)",
    subtitle = "From Theory",
    y = "Theoretical Probability"
  )

Explore Binomial Graphs

• The following arrangement of graphs shows how the shape of binomial distributions changes for \(p \in \{0.1, 0.5, 0.9\}\) as \(n\) increase from 1 to 1024 by powers of 2.
• Notice how the skewness and scale change for each \(p\) as \(n\) increases.

This output is too large to look nice in RStudio, it is formatted nicely in the knitted file.

• When \(p\) is close to the extremes of 0 or 1, the distributions are skewed for small \(n\), but become more symmetric as \(n\) increases.
• When \(p=0.5\), the distributions are always symmetric from the start.
• As the distributions get larger in scale, the same shape emerges: the “normal”, bell-shaped curved.

gbinom(n = 2, p = 0.4)

gbinom(n = 2, p = 0.4) +
  geom_vline(xintercept = 2*0.4, linetype = "dashed", color = "red")

# a and b represent the left and right endpoints of the range you want to graph!
ggplot() +
  geom_binom_density(n = 100, p = 0.5, a = 40, b = 55) +
  geom_binom_density(n = 100, p = 0.6, a = 50, b = 70, color = "red", alpha = 0.5)

Cumulative Probabilities and Quantiles

• In addition to means, variances, and specific probabilities, statistical inference also concerns itself with cumulative probabilities and quantiles.

• Fortunately, R includes quick code to calculate these quantities as well.

• We have already seen dbinom(x, size, prob), which gives \(P(Binom(size, prob) = x)\), and rbinom(n, size, prob), which gives \(n\) random points from the \(Binom(size, prob)\) distribution.

• If you pull up the help page for either function, you will be brought to the general “Binomial Distribution” help page, which brings up four commands. These correspond to the four most common operations for a known probability distribution.

  • These same d<dist>, p<dist>, q<dist>, and r<dist> functions are defined the same way for many distribution types in R!

  • Try ?dnorm for normal distributions, ?dt for “t” distributions, and ?dunif for uniform distributions!

  • There are many other distributions we will not touch on in this class, upon which R supports this core set of four commands, including:

    • F distributions, ?df
    • Chi-squared distributions, ?dchisq
    • Poisson distributions, ?dpois
    • Exponential distributions, ?dexp
    • Gamma and Beta distributions, ?dgamma and ?dbeta

# So we don't have to keep typing out 10 and 0.5
n = 10
p = 0.5

# What is the probability of getting four or fewer heads from 10 coin flips?
pbinom(4, n, p)

## [1] 0.3769531

# The manual alternative:
dbinom(0, n, p) + dbinom(1, n, p) + dbinom(2, n, p) + dbinom(3, n, p) + dbinom(4, n, p)

## [1] 0.3769531

# A shorter, but still manual alternative:
sum(dbinom(0:4, n, p))

## [1] 0.3769531

gbinom(n, p) +
  geom_binom_density(n, p, b = 4, color = "red") +
  scale_x_continuous(breaks = 0:10)

dbinom(4, n, p)

## [1] 0.2050781

# This is P(X > 4), or P(X >= 5)
1 - pbinom(4, n, p)

## [1] 0.6230469

gbinom(n, p) +
  geom_binom_density(n, p, a = 5, color = "red") +
  scale_x_continuous(breaks = 0:10)

# This is P(X > 4), or P(X >= 5). NOT P(X >= 4)!
pbinom(4, n, p, lower.tail = FALSE)

## [1] 0.6230469

# One option: P(X >= 4) is P(X > 3) for the binomial distribution, since it only can be an integer
pbinom(3, n, p, lower.tail = FALSE)

## [1] 0.828125

# Another equally valid option: P(X >= 4) is P(X > 4) + P(X = 4)
pbinom(4, n, p, lower.tail = FALSE) + dbinom(4, n, p)

## [1] 0.828125

gbinom(n, p) +
  geom_binom_density(n, p, a = 4, b = 7, color = "red") +
  scale_x_continuous(breaks = 0:10)

gbinom(n, p) +
  geom_binom_density(n, p, b = 7, color = "red") +
  geom_binom_density(n, p, b = 3, color = "green", alpha = 0.5) +
  scale_x_continuous(breaks = 0:10)

n = 100
p = 0.5
pbinom(56, n, p) - pbinom(51 - 1, n, p)

## [1] 0.3635314

# I use a = 45 to "zoom in" on the calculated area
gbinom(n, p, a = 45, b = 65) +
  geom_binom_density(n, p, a = 45, b = 56, color = "red") +
  geom_binom_density(n, p, a = 45, b = 50, color = "green", alpha = 0.5)

n = 10
p = 0.5

# How much probability is at or below 5 on Binom(10, 0.5)?
pbinom(5, n, p)

## [1] 0.6230469

# How much probability is at or below 6?
pbinom(6, n, p)

## [1] 0.828125

# To capture at least 70% of the probability in Binom(10, 0.5), how large do we need to make x?
qbinom(0.7, n, p)

## [1] 6

# How "far in" do we need to go to capture at least 95% of the probability?
qbinom(0.95, n, p)

## [1] 8

# To capture at least 35% of the probability in Binom(10, 0.5), we need to capture all the bars up to and including...
qbinom(0.35, size = 10, prob = 0.5)

## [1] 4

# Therefore, 4 is the 0.35 quantile of Binom(10, 0.5).

The d, p, q Summary Table

• Recall that dbinom(), pbinom(), and qbinom() are just the binomial example of the larger d, p, q family of commands.

• All of them take a numeric input of interest (as well as the parameters of the distribution; for binom, these are size and prob) and give a numeric output.

• Instead of binom, consider any distribution abbreviation like norm or unif; I’ll write this as <dist>.

• For any <dist>, including binom, these three commands have the following input and output types:

# We continue to use our usual example n and p values, 10 and 0.5.
n = 10
p = 0.5

#dbinom(11, n, p)

#dbinom(0, n, p)

#pbinom(10, n, p)

#pbinom(0, n, p)

#qbinom(1, n, p)

#qbinom(0.5, n, p)

Problem 1

What are the mean and standard deviation of the \(\text{Binomial}(90, 0.7)\) distribution?

# Live coding

Problem 2

Plot the \(\text{Binomial}(90, 0.7)\) distribution. Add a red partly transparent solid vertical line at the mean, dashed vertical lines one standard deviation above and below the mean, and dotted lines two standard deviations above and below the mean.

# Live coding

Problem 3

If \(X \sim \text{Binomial}(90, 0.7)\), what is \(P(X = \mu)\)?

# Live coding

Problem 4

If \(X \sim \text{Binomial}(90, 0.7)\), what is \(P(\mu - \sigma \leq X \leq \mu + \sigma)\)?

# Live Coding

Problem 5

Find the 0.05 and 0.95 quantiles of the \(\text{Binomial}(90, 0.7)\) distribution.

# Live coding

Solutions

n = 90
p = 0.7
mu = n*p
sigma = sqrt(n*p*(1-p))

gbinom(n, p, scale = TRUE) +
  geom_vline(xintercept = mu, color = "red", alpha = 0.5) +
  geom_vline(xintercept = mu + c(-1,1)*sigma,
             color = "red", linetype = "dashed") +
  geom_vline(xintercept = mu + c(-2,2)*sigma,
             color = "red", linetype = "dotted") +
  theme_minimal()

dbinom(mu, n, p)

## [1] 0.09144637

## endpoints
mu + c(-1,1)*sigma

## [1] 58.65259 67.34741

## Within one standard deviation
## P(mu - sigma <= X <= mu + sigma)
## P(59 <= X <= 67) = P(X <= 67) - P(X <= 58)
pbinom(67, n, p) - pbinom(58, n, p)

## [1] 0.6995895

## Sum of probabilities
sum( dbinom(59:67, n, p) )

## [1] 0.6995895

qbinom(c(0.05, 0.95), n, p)

## [1] 56 70

## Compare to the tail probabilities
pbinom(56, n, p) # P(X <= 56) >= 0.05

## [1] 0.06954663

pbinom(55, n, p) # P(X <= 55) < 0.05

## [1] 0.04457106

## Compare to the tail probabilities
pbinom(70, n, p) # P(X <= 70) >= 0.95

## [1] 0.9609796

pbinom(69, n, p) # P(X <= 69) < 0.95

## [1] 0.9354706

The Density Function, Mean, and Variance

• Technically, we should not write \(P(X=x)\) for a continuous variable, because the probability of getting exactly \(x\) is 0, as there are infinitely many possible values. However, we often abuse notation and do anyway.

• Here, we write \(f(x)\), the probability density function. This density function allows for interpretation of probability related to the area under the curve.

\[ f(x) = \frac{1}{\sigma\sqrt{2\pi}} e^{-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2} \]

Explore Normal Graphs

gnorm and geom_norm_density are conceptually identical to gbinom and geom_binom_density above, but require the parameters mu and sigma instead of size and prob.

# Contrasting three normal distributions
ggplot() +
  geom_norm_density(mu = 0, sigma = 1, color = "blue") +
  geom_norm_density(mu = 0, sigma = 0.5, color = "red") +
  geom_norm_density(mu = 2, sigma = 1, color = "darkgreen") +
  labs(
    title = "Three Normal Distributions",
    subtitle = "Blue: N(0, 1). Red: N(0, 0.5). Green: N(2, 1)",
    y = "Density"
  ) +
  geom_hline(yintercept = 0)

Compare the blue to the green curve above. Changing \(\mu\) shifts a normal distribution left or right.

Compare the blue to the red curve above. Changing \(\sigma\) spreads out or narrows a normal distribution.

• Other important features of the normal distribution:
  • No matter what the parameters are, it is a symmetric, “bell-shaped” curve.
    • Unlike the binomial distribution, which we could “skew” for small n’s and extreme values of p, the normal distribution is always symmetric.
  • Normal distributions can exist anywhere along the real number line; positives, negatives, 0, decimal values, et cetera.
    • Technically, all normal distributions have at least some non-zero probability across all values. However, this probability drops to almost 0 once you get more than a few SD’s away from the mean, so we don’t even bother showing it.
  • All normal distributions have total area of 1.0 underneath. That is why the narrower (lower \(\sigma\)) red curve has to be taller than the wider (higher \(\sigma\)) blue and green curves.

geom_norm_fill() instead of geom_norm_density() shades in areas instead of drawing the density line. Optional parameters a and b allow you to specify the left and right endpoints of desired shaded area.

ggplot() +
  geom_norm_fill(mu = 0, sigma = 1, fill = "skyblue") +
  geom_norm_fill(mu = 0, sigma = 1, fill = "blue", alpha = 0.5, a = -0.5, b = 1.5)

dnorm

dnorm(0.5, 0, 1)

## [1] 0.3520653

val = dnorm(0.5, 0, 1)

# An illustration of what dnorm is calculating

# Note: Our custom gnorm() calls the first two arguments "mu" and "sigma", while the d, p, q functions call them "mean" and "sd", but we usually won't name these arguments anyway
gnorm(0, 1, a = -3, b = 3) +
  
  scale_x_continuous(
    breaks = c(-2, 0, 0.5, 2)
  ) +
  
  scale_y_continuous(
    breaks = c(0, 0.2, val, 0.4),
    labels = c("0", "0.2", "dnorm(0.5, 0, 1)", "0.4")
    ) +
  theme(
    panel.grid.minor = element_blank()
  ) +
  
  annotate("segment", x = 0.5, xend = 0.5, y = 0, yend = val, linetype = "solid", size = 2) +
  
  annotate("segment", x = -3, xend = 0.5, y = val, yend = val, linetype = "solid", size = 2)

pnorm()

pnorm(q, mean, sd) takes in a value of the random variable, q, and returns \(P(X <= q)\).

pnorm(1, mean = 0, sd = 1) - pnorm(-0.5, mean = 0, sd = 1)

## [1] 0.5328072

gnorm(0, 1) +
  geom_norm_fill(0, 1, b = 1, fill = "red") +
  geom_norm_fill(0, 1, b = -0.5, fill = "darkgreen", alpha = 0.5)

1 - pnorm(0.7, mean = 0, sd = 1)

## [1] 0.2419637

# or, equivalently:
pnorm(0.7, mean = 0, sd = 1, lower.tail = FALSE)

## [1] 0.2419637

gnorm(100, 25) +
  geom_norm_fill(100, 25, a = NULL, b = 70) +
  geom_norm_fill(100, 25, a = 130, b = NULL) +
  theme_minimal()

left_tail = pnorm(70, 100, 25)

## using subtraction
right_tail = (1 - pnorm(130, 100, 25))
left_tail  + right_tail

## [1] 0.2301393

## using lower.tail
right_tail = pnorm(130, 100, 25, lower.tail = FALSE)
left_tail + right_tail

## [1] 0.2301393

left_tail

## [1] 0.1150697

right_tail

## [1] 0.1150697

# Double the area to the left of 70
2 * pnorm(70, 100, 25)

## [1] 0.2301393

qnorm()

qnorm(p, mean, sd) tells you how large q needs to be for \(P(N(mean, sd) \leq q)\) to be at least p. It is the inverse of pnorm(q, mean, sd).

• How large do we need to make q so \(P(N(100, 25) < q)\) is 70%?

qnorm(0.7, 100, 25)

## [1] 113.11

# I'm not purposefully going to give you the visual this time, I encourage you to try to picture this in your head without it. 
# The N(100, 25) distribution is centered at 100, and we want to capture more than half the area starting from the left. Because the normal distribution is symmetric, we have to go past the center all the way to 113.11 to capture more area.

• How large do we need to make q to capture 15% of the area under \(P(N(100, 25))\)?

qnorm(0.15, 100, 25)

## [1] 74.08917

• Recall the binomial distribution could give you the same value for two different quantiles; for example, 6 is both the 0.7 and the 0.71 quantile of \(Binom(10, 0.5)\).

qbinom(0.7, 10, 0.5)

## [1] 6

qbinom(0.71, 10, 0.5)

## [1] 6

• This will not happen for normal distributions because they are continuous. There is probability (area under the density curve) being added at every single point.

qnorm(0.7, 100, 0.25)

## [1] 100.1311

qnorm(0.71, 100, 0.25)

## [1] 100.1383

mu = 100
sigma = 25

# Answer here
qnorm(0.975, mu, sigma)

## [1] 148.9991

# Answer here
pnorm(75, mu, sigma)

## [1] 0.1586553

# Answer here
2 * pnorm(75, mu, sigma)

## [1] 0.3173105

Standard Normal Distribution and Standardization

The standard normal distribution, often denoted with the random variable \(Z\), is the normal distribution which has mean = 0 and standard deviation = 1; or \(N(0, 1)\).

\(Z\) is a very important statistical distribution whose properties we will study and use significantly.

mu = 100
sigma = 25
x = 80
z = (x-mu)/25

## Using X directly
pnorm(x, mu, sigma)

## [1] 0.2118554

## Standardizing first
## default values for the mean and sd are 0 and 1
pnorm(z)

## [1] 0.2118554

gnorm(mu, sigma) +
  geom_norm_fill(mu, sigma, b=x) +
  theme_minimal()

gnorm() +
  geom_norm_fill(b=z) +
  theme_minimal()

Central Limit Theorem

mu = 100
sigma = 25
gnorm(mu, sigma)

random_data = rnorm(10, mu, sigma)
random_data

##  [1]  51.22708  48.93036  93.15197  80.98691  91.38845  99.54635 109.54453
##  [8]  72.72211 142.43267 102.86998

mean(random_data)

## [1] 89.28004

random_means = tibble(x = replicate(n = 10000, mean(rnorm(10, mu, sigma))))
head(random_means)

## # A tibble: 6 × 1
##       x
##   <dbl>
## 1  85.1
## 2  87.8
## 3  99.4
## 4 100. 
## 5  94.1
## 6  98.5

gnorm(mu, sigma) +
  geom_density(aes(x), random_means) +
  labs(
    title = "Initial (Blue) vs. Sampling Distribution (Black)",
    subtitle = "Sampling Distribution of Sample Mean, n = 10",
    caption = "Notice: Same Center, Less Variance for Sampling Distribution"
  )

# Consistent with theory!
ggplot(random_means, aes(x)) +
  geom_density() +
  geom_norm_density(100, 25/sqrt(10), color = "darkgreen") +
  labs(
    title = "Simulated Dist. of Sample Mean (Black) vs. Theoretical Dist. of Sample Mean (Green)",
    subtitle = "Theoretical Distribution by Central Limit Theorem"
  )

# Consistent with theory!
random_means %>% 
  summarize(mu_simulated = mean(x), sigma_simulated = sd(x))

## # A tibble: 1 × 2
##   mu_simulated sigma_simulated
##          <dbl>           <dbl>
## 1         100.            7.91

random_means_n50 = tibble(x = replicate(n = 10000, mean(rnorm(50, mu, sigma))))

gnorm(mu, sigma) +
  geom_density(aes(x), random_means) +
  geom_density(aes(x), random_means_n50, color = "red") +
  labs(
    title = "Increasing n, Decreasing Variance",
    subtitle = "Initial (Blue) vs. n = 10 Sampling Dist. (Black) vs. n = 50 Sampling Dist. (Red)",
    caption = "Notice: Same Center, Less Variance for Sampling Distribution"
  )

# Consistent with theory!
random_means_n50 %>% 
  summarize(mu_simulated = mean(x), sigma_simulated = sd(x))

## # A tibble: 1 × 2
##   mu_simulated sigma_simulated
##          <dbl>           <dbl>
## 1         99.9            3.51

a = 2
b = 5
random_means_uniform = tibble(x = replicate(n = 10000, mean(runif(10, a, b))))

ggplot(tibble(x=c(a,b), y=1/(b-a))) +
  geom_line(aes(x,y), size=3, color = "blue") +
  geom_density(aes(x), data = random_means_uniform) +
  geom_segment(aes(x, y, xend=x, yend=rep(0,2)), size=3, color = "blue") +
  geom_segment(aes(x, y, x=a-1, y=0, xend=a, yend=0), size=3, color = "blue") +
  geom_segment(aes(x=b+1, y=0, xend=b, yend=0), size=3, color = "blue") +
  scale_x_continuous(breaks = seq(a-1,b+1)) +
  labs(
    title = "Sampling Distribution (Black) of Uniform(2, 5) (Blue)"
  )

# Consistent with theory!
random_means_uniform %>% 
  summarize(mu_simulated = mean(x), sd_simulated = sd(x))

## # A tibble: 1 × 2
##   mu_simulated sd_simulated
##          <dbl>        <dbl>
## 1         3.50        0.270

shape = 0.9
scale = 0.09
exercise_plot = tibble(x = seq(0.01, 1, by = 0.05), density = dgamma(x, shape, scale = scale)) %>% 
  ggplot(aes(x, density)) +
  geom_line()

exercise_plot

# Fill in ???

# exercise_plot +
#   geom_norm_density(mean = ???, sd = ???)

n = 20

xbar_mean = shape * scale
xbar_sd = scale * sqrt(shape) / sqrt(n)


exercise_plot +
  geom_norm_density(mu = xbar_mean, sigma = xbar_sd)

# Your solution here

pnorm(0.1, xbar_mean, xbar_sd)

## [1] 0.8401774

Normal Approximation of Binomial Distributions

• This section is conceptually separate from Central Limit Theorem and only relies on your knowledge of normal and binomial distributions. It will pop up later in statistical inference.

A useful mathematical property of the binomial distribution is that its probability distribution approximately converges to (meaning, “looks like”) a normal distribution if \(n\) is not small and \(p\) is not extremely close to 0 or 1.

• We could see this visually in the “Explore Binomial Graphs” section; as \(n\) got larger for reasonable \(p\) values, the binomial distribution looked more and more like a bell curve.

• Mathematically:

If \(X \sim Binom(n, p)\), then \(X \stackrel{\text{approx}}{\sim}\text{N}(np, \sqrt{np(1-p)})\).

• It is important to remember the binomial distribution is discrete, and the normal distribution is continuous. However, aside from this key difference, their probability distributions are mathematically linked.

• If \(n\) is large enough, and \(p\) is not too close to 0 or 1, the binomial distribution can be approximated by the normal distribution.

• A rule of thumb for “large enough” is if \(np(1-p) \geq 10\).

• Here is an example where the assumption is satisfied.

  • \(n = 100\)
  • \(p = 0.5\)
  • \(np(1-p) = 25\)

## Assumptions satisfied
n = 100
p = 0.5
mu = n*p
sigma = sqrt(n*p*(1-p))

gbinom(n, p, scale=TRUE) +
  geom_norm_density(mu, sigma, color = "red") +
  theme_minimal()

• Here is an example where the assumption is not satisfied.
  • \(n = 100\)
  • \(p = 0.01\)
  • \(np(1-p) = 0.99\)

## Criteria not satisfied
n = 100
p = 0.01
mu = n*p
sigma = sqrt(n*p*(1-p))

gbinom(n, p, scale=TRUE) +
  geom_norm_density(mu, sigma, color = "red") +
  theme_minimal()